Methods and apparatus for automatic magnetic compensation

ABSTRACT

A method for characterizing distortions in the earth&#39;s magnetic field caused by a vehicle having a magnetometer affixed therein is described. The method includes repeatedly measuring the distorted magnetic field utilizing the magnetometer and obtaining a three-dimensional orientation of the vehicle axes with respect to the earth at a time of each magnetometer measurement. The method also includes receiving undistorted earth magnetic field data for the vicinity of the vehicle relative to the earth at the time of each magnetometer measurement and characterizing distortions caused by one or more of the vehicle and magnetometer errors utilizing the magnetic field measurements, the orientations of the vehicle, and the undistorted earth magnetic field data.

CROSS REFERENCE TO RELATED APPLICATIONS

This application claims priority of Provisional Application Ser. No.60/436,980 filed Dec. 30, 2002.

BACKGROUND OF THE INVENTION

This invention relates generally to determination of direction throughmagnetic direction indication, and more specifically to methods andapparatus for compensation of magnetic direction indications to accountfor local magnetic disturbances.

Magnetic direction indicators (e.g., magnetic compasses, magnetometers)are typically compensated to account for local disturbances in anambient magnetic field caused by nearby magnetic objects, for example,ferrous materials, and magnetic fields generated by electrical currents.Known compensation methods are both time consuming and expensiveprocesses. These known compensation methods also provide a limitedaccuracy because the methods are typically optimized for only a singleangle of magnetic inclination.

Some known modem electronic compasses utilize software or firmwarealgorithms for magnetic compensation, but even these devices generallyrequire time-consuming manual compensation procedures to determine a setof optimum compensation coefficients for utilization when performing thealgorithms. Some low precision magnetometers are also known to exist,which provide quick automatic compensation techniques, but thesetechniques only provide a low accuracy compensation. Other knownelectronic compasses utilize biasing circuits as part of a closed loopsystem to attempt to reduce effects of magnetic field disturbances.However, these compasses also incorporate initialization modes, whichcan be complex, and which must be repeated upon each usage of thecompass.

BRIEF SUMMARY OF THE INVENTION

In one aspect, a method for characterizing distortions in the earth'smagnetic field caused by a vehicle is provided. A magnetometer isaffixed to the vehicle and the method comprises repeatedly measuring thedistorted magnetic field utilizing the magnetometer and obtaining athree-dimensional orientation of the vehicle axes with respect to theearth at a time of each magnetometer measurement. The method alsocomprises receiving undistorted earth magnetic field data for thevicinity of the vehicle relative to the earth at the time of eachmagnetometer measurement and characterizing distortions caused by one ormore of the vehicle and magnetometer errors utilizing the magnetic fieldmeasurements, the orientations of the vehicle, and the undistorted earthmagnetic field data.

In another aspect, a method of compensating a magnetometer affixed to avehicle to obtain accurate magnetic heading information for a vehicleorientation is provided. The method comprises using the magnetometer tomeasure a distorted earth magnetic field relative to axes of thevehicle, determining a pitch and roll orientation of the vehicle axeswith respect to the earth, and calculating the distortion of the earth'smagnetic field for any relative angle between the vehicle axes and theearth's undistorted magnetic field. The method further comprisesdetermining a magnetic heading based on the magnetometer measurement,adjusted by the pitch and roll orientation of the vehicle, andcompensated for distortions of the earth's magnetic field.

In yet another aspect, a method for determining a true earth magneticfield from a magnetic field measured by a magnetometer is provided. Themethod comprises generating a truth reference field vector, {tilde over(h)}_(i), from inertial data and a three dimensional map of the earth'smagnetic field, determining a difference between a vector as measured bythe magnetometer and the truth reference vector, and utilizing thedifference to estimate corrections to magnetometer model coefficients.

In still another aspect, a magnetic compass compensation unit fordetermining an orientation of a magnetometer within a vehicle relativeto an undisturbed magnetic field of the earth is provided. Thecompensation unit comprises a processor configured to receivemeasurements of the distorted magnetic field from a magnetometer,receive an orientation of the vehicle axes with respect to the earth ata time corresponding to each magnetometer measurement, receiveundistorted earth magnetic field data for the vicinity of the vehiclerelative to the earth at the time corresponding to each magnetometermeasurement, and characterize the distortions utilizing the magneticfield measurements, the orientations of the vehicle, and the undistortedearth magnetic field data.

In still yet another aspect, a processor programmed to generate a truthreference field vector, {tilde over (h)}_(i), from inertial data and athree dimensional map of the earth's magnetic field is provided. Theprocessor also determines a difference between a vector as measured bythe magnetometer and the truth reference vector, and utilizes thedifference to estimate corrections to magnetometer model coefficients.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is an illustration of a magnetometer in a vehicle passingthrough, and disturbing, a magnetic field.

FIG. 2 is a detailed diagram of a magnetometer.

FIG. 3 is a flow diagram of a method for characterizing distortions inthe earth's magnetic field caused by a vehicle.

FIG. 4 is a flow diagram of a method for compensating a magnetometer forlocal disturbances in the earth's magnetic field.

FIG. 5 depicts an overall structure of a magnetometer autocalibrationand aiding process which utilizes Kalman filtering.

FIG. 6 illustrates magnetometer processing as three functions: ameasurement residual calculation, a measurement matrix calculation, anda measurement covariance matrix calculation.

FIG. 7 details the measurement residual calculation of FIG. 6.

FIG. 8 details the measurement matrix calculation of FIG. 6.

FIG. 9 illustrates a magnetometer process model.

DETAILED DESCRIPTION OF THE INVENTION

Methods and apparatus are herein described which automaticallycompensate a magnetic compass (e.g. a magnetometer) for localdisturbances in the earth's magnetic field caused by a vehicle duringnormal vehicle operations. The methods and apparatus eliminate, orgreatly reduce, time and expense associated with known methods fordetermination of magnetic compensation coefficients. Further, themethods and apparatus provide highly accurate, three dimensionalcompensation coefficients that are valid over all combinations ofvehicle pitch, vehicle roll, vehicle heading, magnetic inclination, andmagnetic declination.

FIG. 1 illustrates an vehicle 10 equipped with a magnetometer 12providing magnetic field strength signals 14. In one embodiment,magnetometer 12 provides magnetic field strength signals 14 representinga measurement of the earth's magnetic field to a magnetic compasscompensation unit 16, which may be part of an aircraft navigationsystem, for example. Magnetometer 12 includes three magnetic sensors(shown in FIG. 2) placed in different orientations from one another. Ina specific embodiment, magnetic sensors are placed such that each sensoris linear with one of the three orthogonal axes of vehicle 10. Magneticfield lines 20 are shown a distance from vehicle 10, and are notdisturbed by presence of vehicle 10. Magnetic field lines 22 however,are disturbed by the presence of vehicle 10. Magnetic field lines 20 and22 are used to illustrate a presence of the earth's magnetic field. Thedisturbance to the earth's magnetic field is illustrated by thenon-uniformities in magnetic field lines 22. The disturbed magneticfield passes through magnetometer 12 which causes sensors to react,based on a magnetic field strength passing across each sensor.Therefore, magnetic field strength signals 14 are inaccurate because ofthe presence of vehicle 10. In one embodiment, magnetometer 12 is asolid state magnetometer.

As is known, magnetic field lines 20 and 22 are representative ofmagnetic field strength. Lines closer together represent a strongermagnetic field strength, and lines farther apart represent a weakermagnetic field. The Earth's magnetic fields originate from the earth'smagnetic north and south poles, and for the most part are uniform.However, and as described above, ferrous materials or other magneticfield sources may interfere with or disturb the earth's magnetic field.Such a disturbance might well be represented by magnetic field lines 22changing direction in the immediate vicinity of a disturbing source,i.e. vehicle 10. Trying to determine magnetic heading of vehicle 10based upon magnetic fields disturbed by vehicle 10 is prone to error.

FIG. 2 is a detailed diagram of magnetometer 12. Magnetometer 12includes an x-axis magnetic sensor 30, a y-axis magnetic sensor 32, anda z-axis magnetic sensor 34 each in a respective orthogonal axis ofvehicle 10. Sensors 30, 32, and 34 generate a signal based on strengthof a magnetic field through which sensors 30, 32, and 34 pass. Sensors30, 32, and 34 are situated orthogonally, and the induced signal will bedifferent for each axes, allowing determination of a three dimensionalorientation of the disturbed magnetic field relative to orientation ofvehicle 10. In one embodiment, magnetic sensors 30, 32, and 34 providesignals to an interface unit 36 which modifies signals from sensors 30,32, and 34 for transmission to magnetic compensation unit 16 (shown inFIG. 1). While sensors 30, 32, and 34 are shown and described as beingorthogonal to one another, the methods and apparatus described hereinare applicable to sensor configurations where at least three sensors areeach oriented in a unique direction, not necessarily orthogonal and notnecessarily along vehicle axes, which allow determination ofcompensation coefficients in a moving vehicle as further describedbelow. Other embodiments exist where only two magnetic sensors areutilized.

FIG. 3 is a flow diagram 50 which illustrates a method performed bymagnetic compass compensation unit 16 to characterize distortions in theearth's magnetic field caused by vehicle 10. The method is utilized todetermine an orientation of magnetometer 12 relative to an undisturbedmagnetic field of the earth, where magnetometer 12 within a vehicle 10which causes a disturbance to the magnetic field of the earth. Referringto flow diagram 50, the distorted magnetic field is repeatedly measured52 utilizing magnetometer 12. A three-dimensional orientation of thevehicle axes is then obtained 54 with respect to the earth at a time ofeach magnetometer measurement. Undistorted earth magnetic field data forthe vicinity of vehicle 10 relative to the earth at the time of eachmagnetometer measurement is received 56. Distortions caused by either orboth of vehicle 10 and magnetometer errors are characterized 58utilizing the magnetic field measurements, the orientations of thevehicle, and the undistorted earth magnetic field data.

In one embodiment, signals 14 from magnetometer 12 are calibrated byunit 16 according to h_(earth)=L_(e) * h_(meas)+h_(pe) where, h_(earth)is a three-dimensional vector for the undistorted Earth's magneticfield, h_(meas) is a three-dimensional vector for the distorted Earth'smagnetic field as measured by magnetometer 12, L_(e) is a three by three(3×3) matrix of magnetic correction coefficients, which includes, forexample, nine correction coefficients, and h_(pe) is a three-dimensionalvector of magnetic correction coefficients which includes, for example,three correction coefficients. Signals 14 from magnetometer 12,h_(meas), are multiplied using a multiplication function by magneticcorrection coefficients, referred to herein as L_(e) which is theaforementioned three by three matrix of magnetic correctioncoefficients. The product of the multiplication is added, using anaddition function, to h_(pe), which is a three-dimensional vector ofmagnetic correction coefficients to provide a magnetic heading which iscompensated for local disturbances to the magnetic field of the earth,for example, caused by the presence of a vehicle.

Magnetic correction coefficients, or L_(e) and h_(pe), are estimated bycomparing the three-dimensional magnetic field strength matrix,h_(meas), as measured by magnetometer 12, against a three-dimensionalmagnetic truth system (not shown). A magnetic truth system provides anindependent means to determine orientation of the Earth's undisturbedmagnetic field relative to the orientation of magnetometer 12. Theundisturbed Earth's magnetic field for a location can be determined froma geographic map of magnetic inclinations and declinations, assumingvehicle position is known through GPS or other means. Pitch and rollorientations of vehicle 10 can be determined by an attitude headingreference system (AHRS). Other sources of pitch and roll signalsinclude, but are not limited to, an attitude reference system, aninertial reference system, and an inertial reference unit. The trueheading orientation of vehicle 10 is provided through utilization of oneof at least four different methods, and then transforming undistortedmagnetic field components to distorted magnetic field components overrelative orientations between the undistorted magnetic field and one ormore of the axes of vehicle 10 and the axes of magnetometer 12.

One method for determining true heading involves a vehicle equipped withGPS. GPS is utilized to determine a true heading while the vehicle ismoving by assuming true heading equals GPS track angle. The methodassumes the vehicle's wheels are aligned with a longitudinal axis of thevehicle and that the vehicle is not turning or side slipping. However,effects of turns during movement are easily compensated by knowing thegeometry of, for example, aircraft landing gears and by measuring GPSvelocities and/or the yaw rate as determined by aircraft gyroscopes.

A second method for determining true heading involves a vehicle equippedwith GPS and an Attitude and Heading Reference System (AHRS). Well-knownAHRS algorithms are utilized to calculate a true heading while thevehicle is experiencing horizontal accelerations, including centrifugalaccelerations caused by turns. These AHRS algorithms can continue todetermine a true heading for a limited time duration after theacceleration stops, on the order of minutes.

If the vehicle has an Inertial Reference System (IRS), then a thirdmethod for determining true heading involves use of a true headingsignal received directly from the IRS. An alternative is a synthesizedmagnetic heading signal from the IRS. A fourth method for determiningtrue heading is utilization of a dual antenna GPS, from which trueheading can be directly calculated.

Twelve coefficients, including nine magnetic correction coefficients andthree vector of magnetic correction coefficients, are calculated usingseveral methods. A preferred embodiment is to populate athree-dimensional magnetic error table with data during normal aircraftoperations. The magnetic error table, in one embodiment, maintainsseparate data points for multiple orientations of the earth's actualmagnetic field with respect to vehicle pitch, roll, and yaw axes. Forexample, the magnetic error table may include separate data points at 30degree increments in yaw and at 60 degree increments in pitch and roll.In such an embodiment, the magnetic error table is configured with 216data points. Each data point is a three-dimensional vector of magneticerror relative to the vehicle axes. Therefore, magnetic error tableincludes 648 scalar values. The data at each data point, in oneembodiment, is low pass filtered, with a time constant of several hours.The low pass filtering occurs during a period when the orientation ofthe actual magnetic field relative to the orientation of vehicle 10corresponds to the three-dimensional vector as represented by one ofdata points. Data points are not updated when the magnetic field presentat vehicle 10 does not correspond to a particular data point.

Several, but perhaps not all, of the data points in magnetic error tablewill be updated during each flight, depending on, for example, how manyturns vehicle 10 makes, vehicle orientation, and the angles ofinclination it flies through. After a number of flights, the twelvecompensation coefficients are determined from the highly overdetermined648 scalar values of the 216 data points, through well-knownmathematical procedures, for example, utilizing a least squaresdetermination.

In conjunction with the above, a method for compensating a magnetometeraffixed to a vehicle is illustrated in flowchart 60 of FIG. 4. Themethod provides accurate magnetic heading information for a vehicleorientation. Referring to flowchart 60, a magnetometer 12 (shown inFIG. 1) is used 62 to measure a distorted earth magnetic field relativeto axes of vehicle 10 (shown in FIG. 1), and a pitch and rollorientation of the axes of vehicle 10 is determined 64 with respect tothe earth. A distortion of the earth's magnetic field for any relativeangle between the vehicle axes and the earth's undistorted magneticfield is calculated 66 and a magnetic heading based on the magnetometermeasurement, adjusted by the pitch and roll orientation of vehicle 10,and compensated for distortions of the earth's magnetic field isdetermined 68.

In a further embodiment, a Kalman filter model is developed to carry outthe above described procedures, and which automatically calibrates thetwelve parameters that make up the three dimensional magnetometer modelfor estimating the true earth's magnetic field from the field measuredby magnetometer 12. This calibration method uses pitch, roll, heading,and position data from one or more of a GPS, AHRS, inertial referencesystem, inertial navigation system and ground based navigational aids(i.e., VOR, distance measuring equipment (DME)), along with a threedimensional map of the earth's magnetic field to generate the “truth”reference field vector. The difference between the magnetometer-measuredfield vector and the “truth” reference vector becomes the measurementfor the Kalman filter. The Kalman filter then estimates corrections tothe magnetometer model coefficients as well as corrections to theattitude and heading angles of the navigation system.

FIG. 5 depicts an overall structure 100 of the above describedmagnetometer auto-calibration and compensation process. Magnetometer 102measures the earth's magnetic field and provides the measurements to amagnetometer measurement process 104. A GPS receiver 106 provides GPSmeasurements to a GPS measurement process 108. Although GPS receiver 106is the only aiding source depicted in FIG. 4 besides magnetometer 102,other aiding sources may be utilized as well, such as barometericaltitude, airspeed, and other sources as previously described. Outputsfrom magnetometer measurement process 104 and GPS measurement process108 are combined utilizing measurement model 110, whose outputs areinput to Kalman filter 112.

Kalman filter 112 provides a filter state output 114 which is providedas an input to magnetometer measurement process 104 and GPS measurementprocess 108 as well as to navigation solution generator 116, whichprovides, in the embodiment shown, an aided navigation solution 118.Inertial sensors 120 also provide data to navigation solution generator116. The aided navigation solution 118 is fed back as inputs tomagnetometer measurement process 104 and GPS measurement process 108 aswell as to a GPS/inertial process model 122 and a magnetometer processmodel 124. Outputs of GPS/inertial process model 122 and magnetometerprocess model 124 are input to a process model 126, whose output isinput to Kalman filter 112.

In developing the Kalman filter model, a model of the field as measuredby magnetometer 12 is utilized. Specifically, the field measured by themagnetometer {tilde over (h)}_(m) is related to the true earth's field{tilde over (h)} according to:{tilde over (h)} _(m) =M{tilde over (h)}+ {tilde over (h)} _(p)   (1)

where

-   -   M=3×3 magnetic permeability matrix    -   {tilde over (h)}_(p)=3×1 vector of field offset errors resulting        from permanent magnetization

A tilde ({tilde over ()}) is utilized to distinguish the magnetic fieldvectors from the Kalman filter measurement vector h which is describedbelow.

By solving for {tilde over (h)} the true field from the measured fieldis computed via $\begin{matrix}\begin{matrix}{{\overset{\sim}{h} = {M^{- 1}\left( {{\overset{\sim}{h}}_{m} - {\overset{\sim}{h}}_{p}} \right)}}\quad} \\{= {L\left( {{\overset{\sim}{h}}_{m} - {\overset{\sim}{h}}_{p}} \right)}}\end{matrix} & (2)\end{matrix}$

where the matrix L is defined as L=M⁻¹. Comparing equation (2) abovewith the calibration equation from the first embodiment, (i.e.h_(earth)=L_(e) * h_(meas)+h_(pe)), it is seen that the two equationsare equivalent if L_(e)=L and h_(pe)=−L{tilde over (h)}_(p). To developthe Kalman filter, it is assumed that the initial value of matrix L isequal to the identity matrix, and that the initial vector of fieldoffset errors, {tilde over (h)}_(p), is zero.

A linearized error equation for the magnetometer-measured andcompensated body-axis earth's field components is obtained by takingpartial differentials of equation (2). For example, $\begin{matrix}{{\begin{matrix}{{\delta\overset{\sim}{h}} = {{\delta\quad{L\left( {{\overset{\sim}{h}}_{m} - {\overset{\sim}{h}}_{p}} \right)}} - {L\quad\delta{\overset{\sim}{h}}_{p}}}} \\{= {{\delta\quad L\quad\Delta\quad\overset{\sim}{h}} - {L\quad\delta{\overset{\sim}{h}}_{p}}}}\end{matrix}\quad{{where}\quad\Delta\quad\overset{\sim}{h}\quad{is}\quad{defined}\quad{as}}{{\Delta\quad\overset{\sim}{h}} = {{\overset{\sim}{h}}_{m} - {{\overset{\sim}{h}}_{p}.}}}}\quad} & (3)\end{matrix}$

Equation (3) provides a basic estimate of errors in a matrix format andis broken into three separate measurement equations (one for each fieldcomponent), where each equation is in the form of a row (measurementmapping) vector multiplied by a column vector of magnetometer parametererrors to provide scalars that are functions of the vectors, forexample,δ{tilde over (h)} _(x) =Δ{tilde over (h)} ^(T) δl _(r1) −l _(r1) ^(T)δ{tilde over (h)} _(p)   (4)δ{tilde over (h)} _(y) =Δ{tilde over (h)} ^(T) δl _(r2) −l _(r2) ^(T)δ{tilde over (h)} _(p)   (5)δ{tilde over (h)} _(z) =Δ{tilde over (h)} ^(T) δl _(r3) −l _(r3) ^(T)δ{tilde over (h)} _(p)   (6)

where

-   -   δl_(ri) ^(T)=the i^(th) row of δL    -   l_(ri) ^(T)=the i^(th) row of L

Alternatively, each term is written as a matrix times an error vector bydefining the following 9×1 error vector $\begin{matrix}{{{\delta\quad I} = \begin{bmatrix}{\delta\quad I_{r1}} \\{\delta\quad I_{r2}} \\{\delta\quad I_{r3}}\end{bmatrix}}\quad} & (7)\end{matrix}$

Then, equation (3) can be re-written as $\begin{matrix}{{{\delta\overset{\sim}{h}} = {{\begin{bmatrix}{\Delta\quad{\overset{\sim}{h}}^{T}} & 0_{1 \times 3} & 0_{1 \times 3} \\0_{1 \times 3} & {\Delta\quad{\overset{\sim}{h}}^{T}} & 0 \\0_{1 \times 3} & 0_{1 \times 3} & {\Delta\quad{\overset{\sim}{h}}^{T}}\end{bmatrix}\quad\delta\quad I} - {L\quad\delta\quad{\overset{\sim}{h}}_{p}}}}\quad} & (8)\end{matrix}$

which provides magnetometer errors.

To compute a GPS/AHRS-based or inertial navigation system-based “truth”field {tilde over (h)}_(i), the earth's field vector is first determinedin north/east/down frame components {tilde over (h)}^(N) based on thecurrent latitude, longitude, altitude, and (perhaps) time using anaccurate model or map. Next the earth's field vector is transformed intobody coordinates, to find errors in the “truth” source via the followingtransformations:{tilde over (h)}_(i)=C_(L) ^(B)C_(N) ^(L){tilde over (h)}^(N)   (9)

where C_(L) ^(B) transforms a vector from the local-level frame(L-frame) to the body frame (B-frame) and is the transpose of theattitude direction cosine matrix C, and C_(N) ^(L) accounts for therotation in azimuth of the local-level frame with respect to north bythe wander angle α and is given by ${C_{N}^{L} = \begin{bmatrix}{\cos\quad\alpha} & {\sin\quad\alpha} & 0 \\{{- \sin}\quad\alpha} & {\cos\quad\alpha} & 0 \\0 & 0 & 1\end{bmatrix}}\quad$

By taking partial differentials of equation (9), $\begin{matrix}\begin{matrix}{{\delta{\overset{\sim}{h}}_{i}} = {{\delta\quad C^{T}C_{N}^{L}{\overset{\sim}{h}}^{N}} + {C^{T}\delta\quad C_{N}^{L}{\overset{\sim}{h}}^{N}} + {C^{T}C_{N}^{L}\delta{\overset{\sim}{h}}^{N}}}} \\{{= {{\left\lbrack {C^{T}\left\{ \gamma \right\}} \right\rbrack C_{N}^{L}{\overset{\sim}{h}}^{N}} + {{C^{T}\left\lbrack {{- \left\{ ɛ \right\}}C_{N}^{L}} \right\rbrack}{\overset{\sim}{h}}^{N}} + {C^{T}C_{N}^{L}\delta{\overset{\sim}{h}}^{N}}}}\quad}\end{matrix} & (10)\end{matrix}$

where γ is the attitude error vector which represents the angular errorof the L-frame relative to the B-frame, ε the angular position errorvector which represents the angular error of the L-frame relative to theearth frame (E-frame), and {v} represents the skew-symmetric matrix formof the enclosed vector v and is defined by $\begin{matrix}{{\left\{ v \right\} \equiv \begin{bmatrix}0 & {- v_{z}} & {\quad v_{y}} \\v_{z} & 0 & {- v_{x}} \\{- v_{y}} & v_{x} & 0\end{bmatrix}}\quad} & (11)\end{matrix}$

During the GPS-aided mode, the “psi-angle” inertial error model isimplemented in the Kalman filter, shown below. In such an embodiment,the attitude error states are actually the three components of theangular error vector Ψ defined asΨ=γ−ε  (12)

Using this substitution in equation (10), the first two terms combine togive yield followingδ{tilde over (h)} _(i) =C ^(T) {Ψ}C _(N) ^(L) {tilde over (h)} ^(N) +C^(T) C _(N) ^(L) δ{tilde over (h)} ^(N)   (13)

To express each term in equation (13) as a matrix times an error vector,the first term must be re-arranged. The product of a skew-symmetricmatrix form of a vector v multiplied by a vector u is the same as thecross product v x u. By reversing the order of the cross-product, thesign is reversed. Therefore,{Ψ}C _(N) ^(L) {tilde over (h)} ^(N) =Ψx(C _(N) ^(L){tilde over(h)}^(N))=−(C _(N) ^(L) {tilde over (h)} ^(N))xΨ=−{C _(N) ^(L) {tildeover (h)} ^(N)}Ψ  (14)

and equation (13) is rewritten asδ{tilde over (h)} _(i) =−C ^(T) {C _(N) ^(L) {tilde over (h)} ^(N) }Ψ+C^(T) C _(N) ^(L) δ{tilde over (h)} ^(N)   (15)

The measurements to the Kalman filter are the differences between eachof the magnetometer-derived earth's field components and the respectiveGPS/AHRS-derived earth's field components. These measurements can beprocessed one at a time or all at once. Each measurement equation is inthe form ofz _(k) =h _(k) ^(T) x _(k) +v _(k)   (16)

Where h_(k) is the measurement mapping vector at time t_(k), x_(k) isthe state vector, and v_(k) is white (uncorrelated) measurement noise.If all are processed at once, the measurement equation takes the form ofz _(k) =H _(k) x _(k) +v _(k)   (17)

where z_(k) is a 3×1 vector of measurements at time t_(k), H_(k) is a3×n measurement matrix, and v_(k) is the measurement noise vector.

By assuming that the three “psi-angle” attitude errors make up firstthree states and the twelve magnetometer parameters and the three earthfield modeling errors make up the last fifteen states. The state vectoris then written asx^(T)=[Ψ^(T) . . . other nav/sensor states . . . δ1 ^(T) δ{tilde over(h)}_(p) ^(T) (δ{tilde over (h)}^(N))^(T)]  (18)

The measurement vector is determined subtracting equation (15) fromequation (8) and adding measurement noise $\begin{matrix}\begin{matrix}{z = {{\overset{\sim}{h} - {\overset{\sim}{h}}_{i}} = {{\delta\overset{\sim}{h}} - {\delta{\overset{\sim}{h}}_{i}} + v}}} \\{{= {{\begin{bmatrix}{\Delta\quad{\overset{\sim}{h}}^{T}} & 0_{1 \times 3} & 0_{1 \times 3} \\0_{1 \times 3} & {\Delta\quad{\overset{\sim}{h}}^{T}} & 0 \\0_{1 \times 3} & 0_{1 \times 3} & {\Delta\quad{\overset{\sim}{h}}^{T}}\end{bmatrix}\quad\delta\quad I} - {L\quad\delta\quad{\overset{\sim}{h}}_{p}} + {C^{T}\left\{ {C_{N}^{L}\quad{\overset{\sim}{h}}^{N}} \right\}\psi} - {C^{T}C_{N}^{L}\delta{\overset{\sim}{h}}^{N}} + v}}\quad}\end{matrix} & (19)\end{matrix}$

where v represents a vector of uncorrelated measurement noise due toelectromagnetic noise and random earth field modeling errors such asquantization.

From equations (17), (18), and (19) it is shown that the measurementmapping matrix H is $\begin{matrix}\begin{matrix}{{H = \begin{bmatrix}H_{\psi} & 0 & \cdots & 0 & H_{\delta 1} & H_{{\overset{\sim}{h}}_{p}} & H_{{\overset{\sim}{h}}^{N}}\end{bmatrix}}\quad} \\{where} \\{H_{\psi} = {C^{T}\left\{ {C_{N}^{L}\quad{\overset{\sim}{h}}^{N}} \right\}}} \\{H_{\delta 1} = \begin{bmatrix}{\Delta\quad{\overset{\sim}{h}}^{T}} & 0_{1 \times 3} & 0_{1 \times 3} \\0_{1 \times 3} & {\Delta\quad{\overset{\sim}{h}}^{T}} & 0 \\0_{1 \times 3} & 0_{1 \times 3} & {\Delta\quad{\overset{\sim}{h}}^{T}}\end{bmatrix}} \\{H_{{\overset{\sim}{h}}_{p}} = {- L}} \\{H_{{\overset{\sim}{h}}^{N}} = {{- C^{T}}C_{N}^{L}}}\end{matrix} & (20)\end{matrix}$

Barring any disturbance (such as lightning strikes), the magnetometerparameters are assumed to be very stable over time, therefore areasonable model to describe random processes for each of these statesis a Gauss-Markov process with a very long correlation time, in oneembodiment, about 1000 hours. A reasonable one-sigma value for each ofthe nine δL states is approximately 0.15 (15%), which corresponds to anangular error uncertainty of about 8.5°. A reasonable one-sigma valuefor each of the three offset errors contained in the vector δh_(p) is 75mgauss. Assuming an earth's field of 500 mgauss, this would againcorrespond to an angular error of about 8.5°. (The earth's fieldintensity ranges from about 250-650 mgauss). If a major disturbance doesoccur, it would be rather abrupt and would certainly be detected bystandard Kalman filter residual screening.

A first-order Gauss-Markov process x is defined by the following$\begin{matrix}{{\overset{.}{x} = {{{- \frac{1}{\tau}}x} + {\sqrt{\frac{2\quad\sigma^{2}}{\tau}}{u(t)}}}}\quad} & (21)\end{matrix}$

where u(t) is unity white noise. The discrete-time model of this processis given byx _(k+1)=Φ_(k) x _(k) +w _(k)   (22)

where

-   -   Φ_(k)=e^(−T/) ^(τ)    -   w_(k)=Gaussian white noise sequence with variance=σ²(1−e^(−2T/)        ^(τ) )    -   σ=steady-state one-sigma value of the process    -   τ=correlation time of the process    -   T=the discrete time interval between t_(k) and t_(k+1)

The earth field modeling errors will vary as a function of positionrather than time. Therefore, a Gauss-Markov process model which variesover distance traveled would be appropriate as such a model defines howstate vectors vary over time. According to the National Geophysical DataCenter (NGDC), local anomalies in geomagnetism which deviate from theInternational Geomagnetic Reference Field (IGRF) can exceed 10° ofinclination and/or declination. In fact, Minnesota contains an extremedeclination anomaly of 16° east in one area and 12° west just a fewmiles away—a change of 28 degrees. Over most of the world, however, theIRGF model is good to 0.5°. A reasonable one-sigma value to use to boundthese larger anomalies might be 50 mgauss (corresponding to a 5.7 degreeangular error for a 500 mgauss field) with a correlation distance of 10nautical miles. Thus when cruising at 500 knots, the effectivecorrelation time would be 72 seconds. However, during final approach at150 knots an effective correlation time would be about four minutes.

The one-sigma measurement noise for the magnetometer is assumed to be onthe order of 0.5 mgauss for each component which corresponds to anangular error of 0.06° for a 500 mgauss field.

The GPS/inertial process models are augmented by the fifteen processmodels associated with the magnetometer states to form the overallprocess model in matrix formx _(k+1) =Φ _(k) x _(k) +w _(k)   (23)

where

-   -   Φ_(k)=state transition matrix=the undriven response over time        t_(k) to t_(k+1)    -   w_(k)=process noise vector=driven response to the white noise        input over time t_(k) to t_(k+1)

The magnetometer measurements may be processed independently or alongwith other aiding measurements such as GPS by augmenting the measurementvector z with these three additional measurements. Likewise, themeasurement matrix H would also be augmented with these three additionalrows. Standard Kalman filter time propagation and measurement updateequations are iterated using known methods. For example, the estimatedstate vector {circumflex over (x)}_(k) and its associated error noisecovariance matrix P_(k) are projected forward to the next time step asfollows{circumflex over (x)}_(k+1) ⁻=Φ_(k){circumflex over (x)}_(k)   (24)P _(k+1) ⁻=Φ_(k) P _(k)Φ_(k) ^(T) +Q _(k)   (25)

where

-   -   {circumflex over (x)}_(k+1)=estimated error state vector prior        to the measurement update at time t_(k+1)    -   P_(k)=E[{circumflex over (x)}_(k){circumflex over (x)}_(k)        ^(T)]=error covariance matrix after the measurement update at        time t_(k)    -   Q_(k)=E [w_(k)w_(k) ^(T)]=process noise covariance matrix over        time t_(k) to t_(k+1)    -   P_(k+1) ⁻=the error covariance matrix prior to the measurement        update at time t_(k+1)

and the estimated state vector and the error covariance matrix areupdated for the new measurement vector z_(k) according toK _(k) =P _(k) ⁻ H _(k) ^(T)(H _(k) P _(k) ⁻ H _(k) ^(T) +R _(k))⁻¹  (26){circumflex over (x)} _(k) ={circumflex over (x)} _(k) ⁻ +K _(k)(z _(k)−H _(k) {circumflex over (x)} _(k) ⁻)   (27)P _(k)=(I−K _(k) H _(k))P _(k) ⁻  (28)

where

-   -   R_(k)=E[v_(k)v_(k) ^(T)]=measurement noise covariance matrix at        time t_(k).

FIG. 6 is a detailed illustration of one embodiment of magnetometermeasurement process 104, which includes three functions, a measurementresidual calculation 140, a measurement matrix calculation 142, and ameasurement covariance matrix calculation 144. To calculate measurementresiduals, measurement residual calculation 140 receives as inputs,measured earth's magnetic field as measured by magnetometer 102 (shownin FIG. 5), filter state output 114 from Kalman filter 112 (shown inFIG. 5) and aided navigation solution 118 (also shown in FIG. 5).Measurement matrix calculation 142 utilizes results of residualcalculations made in residual calculation 140 to provide a magneticfield matrix output 146. Measurement covariance matrix calculation 144is based on an identity matrix 148.

FIG. 7 illustrates measurement residual calculation 140 (also shown inFIG. 6). Measured earth magnetic field 160 is first compensated usingthe current estimates of δ{tilde over (h)}_(p) 162 and L 164, an inverseof the permeability matrix M. Note that a nominal value for L 164 isassumed to be identity matrix 166. The matrix δL 168 is an estimate ofthe error from this nominal identity matrix 166 and an estimate of theerrors from magnetometer 12, δL, are subtracted from it to yield acurrent estimate of L 164. The magnetic “truth” field 170 (an inertialmodel of the earth's field) is computed from current estimates ofaltitude 172, position 174, and current date and time 176 using earth'sfield model 178 and current earth-to-local-level frame andlocal-level-to-body frame transformations 180 and is corrected for theearth field modeling error estimate δ{tilde over (h)}^(N) 182. Since themeasured field vector 184 and truth reference vector 186 are correctedexplicitly (and the attitude error is reset in the strapdown navigationequations) with the state vector estimate from the filter, thedifference between the two corrected field vectors 184 and 186 becomesthe measurement residual vector 188 (z_(k)−H_(k){circumflex over(x)}_(k) ⁻) of equation (27).

FIG. 8 illustrates measurement matrix calculation 142 which provides asolution of equation (20), which is described above. Current date andtime 176 is input as are local level transformations 180 and a modeledearth magnetic field 178. Also input to measurement matrix calculation142 are earth field modeling error estimate δ{tilde over (h)}^(N) 182and L 164, an inverse of the permeability matrix M. An output ofmeasurement matrix calculation 142 is the measurement mapping matrix H.

FIG. 9 shows a magnetometer process model 200 which results incalculation of Φ_(k), the state transition matrix described with respectto equation (23) above, and Q_(k) the process noise covariance matrixover time, as described with respect to equation (25), for themagnetometer states.

Once the twelve compensation components are determined, either throughKalman filtering, or through population of the magnetic error table,unit 16 is configured to remove local disturbances to the earth'smagnetic field. Once the local disturbances are removed, accuratenavigation of an aircraft or other vehicle is attained since themagnetic heading information supplied to the pilot is a true magneticheading, uncompromised by local magnetic disturbances.

The above described methods and apparatus illustrate how magnetometerreadings can be automatically compensated to remove local disturbancesin the earth's magnetic field based upon data received from, forexample, an inertial navigation system. While the invention has beendescribed in terms of various specific embodiments, those skilled in theart will recognize that the invention can be practiced with modificationwithin the spirit and scope of the claims.

1. A method for characterizing distortions in the earth's magnetic fieldcaused by a vehicle, a magnetometer affixed to the vehicle, said methodcomprising: repeatedly measuring the distorted magnetic field utilizingthe magnetometer; obtaining a three-dimensional orientation of thevehicle axes with respect to the earth at a time of each magnetometermeasurement; receiving undistorted earth magnetic field data for thevicinity of the vehicle relative to the earth at the time of eachmagnetometer measurement; and utilizing the magnetic field measurements,the orientations of the vehicle, and the undistorted earth magneticfield data to characterize distortions caused by one or more of thevehicle and magnetometer errors.
 2. A method according to claim 1wherein obtaining an orientation of the vehicle axes comprises receivingpitch and roll signals from one or more of an attitude heading referencesystem, an attitude reference system, an inertial reference system, andan inertial reference unit.
 3. A method according to claim 2 whereinobtaining an orientation of the vehicle axes comprises: moving a vehicleequipped with a GPS; determining a GPS tracking angle; and assuming theGPS tracking angle is the true heading.
 4. A method according to claim 2wherein obtaining an orientation of the vehicle axes comprises: moving avehicle equipped with a GPS; and calculating a true heading during, andimmediately following, vehicle turns or accelerations, utilizingGPS-augmented attitude and heading reference system algorithms.
 5. Amethod according to claim 1 wherein obtaining an orientation of thevehicle axes comprises receiving pitch, roll, and at least one of a trueheading signal and a synthetic magnetic heading signal from an inertialreference system.
 6. A method according to claim 1 wherein obtaining anorientation of the vehicle comprises: receiving pitch and roll signalsfrom an attitude heading reference system; and receiving at least one ofa true heading signal and a synthetic magnetic heading signal from aninertial reference system.
 7. A method according to claim 1 whereinobtaining an orientation of the vehicle comprises: receiving pitch androll signals from an attitude heading reference system; and determininga true heading utilizing a dual antenna GPS.
 8. A method according toclaim 1 wherein receiving undistorted earth magnetic field datacomprises determining vehicle position and using a table of magneticinclination and declination data for a surface of the earth.
 9. A methodaccording to claim 1 wherein characterizing distortions caused by thevehicle comprises transforming undistorted magnetic field components todistorted magnetic field components over relative orientations betweenthe undistorted magnetic field and one or more of the axes of thevehicle and the axes of the magnetometer.
 10. A method according toclaim 1 wherein obtaining a three-dimensional orientation of the vehicleaxes comprises determining vehicle position using at least one of aninertial navigational system and ground based navigation aids.
 11. Amethod according to claim 1 wherein characterizing distortions comprisescalculating correction coefficients for the magnetometer.
 12. A methodaccording to claim 11 wherein calculating correction coefficientscomprises compensating the magnetometer by estimating L_(e) and H_(pe)in the function H_(earth)=L_(e) * H_(meas)+H_(pe) where, H_(earth) is athree-dimensional vector representing the undisturbed Earth's magneticfield, H_(meas) represents a three-dimensional vector for the disturbedEarth's magnetic field as measured by magnetometer, L_(e) is a three bythree (3×3) matrix of magnetic correction coefficients, and H_(pe) is athree-dimensional vector of magnetic correction coefficients.
 13. Amethod of compensating a magnetometer affixed to a vehicle to obtainaccurate magnetic heading information for a vehicle orientation, saidmethod comprising using the magnetometer to measure a distorted earthmagnetic field relative to axes of the vehicle; determining a pitch androll orientation of the vehicle axes with respect to the earth;calculating the distortion of the earth's magnetic field for anyrelative angle between the vehicle axes and the earth's undistortedmagnetic field; and determining a magnetic heading based on themagnetometer measurement, adjusted by the pitch and roll orientation ofthe vehicle, and compensated for distortions of the earth's magneticfield.
 14. A method according to claim 13 wherein determining a magneticheading further comprises compensating the magnetometer by solving thefunction H_(earth)=L_(e) * H_(meas)+H_(pe) by estimating L_(e) andH_(pe) where, H_(earth) is a three-dimensional vector representing theundisturbed Earth's magnetic field, H_(meas) represents athree-dimensional vector for the disturbed Earth's magnetic field asmeasured by magnetometer, L_(e) is a three by three (3×3) matrix ofmagnetic correction coefficients, and H_(pe) is a three-dimensionalvector of magnetic correction coefficients.
 15. A method according toclaim 13 further comprising augmenting pitch and roll orientation databased on magnetometer measurements and the magnetometer compensation.16. A method according to claim 13 further comprising estimating L_(e)and H_(pe) based on multiple measurements of H_(earth) and H_(meas). 17.A method according to claim 16 wherein determining a magnetic headingcomprises: determining vehicle position and using a table of magneticinclination and declination data for a surface of the earth; comparingthe measured magnetic field, H_(meas), against the magnetic inclinationand declination data to determine the magnetic correction coefficients,H_(pe); and providing L_(e) and H_(pe) as correction coefficients to themagnetometer.
 18. A method according to claim 13 further comprisingdetermining the vector of magnetic correction coefficients, H_(pe), andthe matrix of magnetic correction coefficients, L_(e), utilizing athree-dimensional magnetic error table, the error table configured tomaintain separate data points for multiple aircraft orientations betweenan actual magnetic field and vehicle pitch, roll, and yaw axes.
 19. Amethod according to claim 18 wherein the magnetic error table comprisesseparate data points at 30 degree increments in yaw and at 60 degreeincrements in pitch and roll, each data point being a three-dimensionalvector of magnetic error relative to vehicle axes.
 20. A methodaccording to claim 19 comprising determining the correction coefficientsutilizing scalar values from the data points.
 21. A method fordetermining a true earth magnetic field from a magnetic field measuredby a magnetometer, said method comprising: generating a truth referencefield vector, {tilde over (h)}_(i), from sources of pitch, roll,heading, and position independent of the magnetometer and a threedimensional map of the earth's magnetic field; determining a differencebetween a vector as measured by the magnetometer and the truth referencevector; and utilizing the difference to estimate corrections tomagnetometer model coefficients.
 22. A method according to claim 21further comprising utilizing the difference to estimate corrections toone or more of vehicle pitch, roll, and heading.
 23. A method accordingto claim 21 wherein pitch, roll, and heading data is provided by atleast one of a GPS, an attitude and heading reference system, aninertial reference system, an inertial navigation system, and groundbased navigational aids.
 24. A method according to claim 21 furthercomprising utilizing the difference to provide corrections to attitudeand heading angles of a navigation system.
 25. A method according toclaim 21 further comprising estimating corrections using a Kalmanfilter.
 26. A method according to claim 21 wherein generating a truthreference field vector comprises relating the field measured by themagnetometer {tilde over (h)}_(m) to the true earth's field {tilde over(h)} according to {tilde over (h)}_(m)=M{tilde over (h)}+{tilde over(h)}_(p), where M is a 3×3 magnetic permeability matrix and {tilde over(h)}_(p) is a 3×1 vector of field offset errors resulting from permanentmagnetization.
 27. A method according to claim 26 wherein determining adifference comprises: solving for the true earth's field, {tilde over(h)}, according to ${\begin{matrix}{\overset{\sim}{h} = {M^{- 1}\left( {{\overset{\sim}{h}}_{m} - {\overset{\sim}{h}}_{p}} \right)}} \\{= {L\left( {{\overset{\sim}{h}}_{m} - {\overset{\sim}{h}}_{p}} \right)}}\end{matrix},}\quad$ where the matrix L is defined as L=M⁻¹; andassuming that the initial value of matrix L is equal to the identitymatrix and the initial vector of field offset errors, {tilde over(h)}_(p), is zero.
 28. A method according to claim 27 furthercomprising: obtaining a linearized error equation for themagnetometer-measured and compensated body-axis earth's field componentsaccording to ${\begin{matrix}{{\delta\overset{\sim}{h}} = {{\delta\quad{L\left( {{\overset{\sim}{h}}_{m} - {\overset{\sim}{h}}_{p}} \right)}} - {L\quad\delta{\overset{\sim}{h}}_{p}}}} \\{= {{\delta\quad L\quad\Delta\quad\overset{\sim}{h}} - {L\quad\delta{\overset{\sim}{h}}_{p}}}}\end{matrix},}\quad$ where Δ{tilde over (h)} is defined as Δ{tilde over(h)}={tilde over (h)}_(m)−{tilde over (h)}_(p); and forming a separatemeasurement equation for each field component from the linearized errorequation according to δ{tilde over (h)}_(x)=Δ{tilde over(h)}^(T)δl_(r1)−l_(r1) ^(T)δ{tilde over (h)}_(p), δ{tilde over(h)}_(y)=Δ{tilde over (h)}^(T)δl_(r2)−l_(r2) ^(T)δ{tilde over (h)}_(p),and δ{tilde over (h)}_(z)=Δ{tilde over (h)}^(T)δl_(r3)−l_(r3)^(T)δ{tilde over (h)}_(p), where each equation is in the form of a rowvector multiplied by a column vector of magnetometer parameter errors toprovide scalars that are functions of the vectors.
 29. A methodaccording to claim 28 where magnetometer errors are provided accordingto ${{\delta\overset{\sim}{h}} = {{\begin{bmatrix}{\Delta\quad{\overset{\sim}{h}}^{T}} & 0_{1 \times 3} & 0_{1 \times 3} \\0_{1 \times 3} & {\Delta\quad{\overset{\sim}{h}}^{T}} & 0 \\0_{1 \times 3} & 0_{1 \times 3} & {\Delta\quad{\overset{\sim}{h}}^{T}}\end{bmatrix}\quad\delta\quad I} - {L\quad\delta\quad{{\overset{\sim}{h}}_{p}.}}}}\quad$30. A method according to claim 21 wherein generating a truth referencefield vector, {tilde over (h)}_(i), from inertial data and a threedimensional map of the earth's magnetic field comprises: determining theearth's field vector in north/east/down frame components {tilde over(h)}^(N) based on one or more of a current latitude, longitude,altitude, and the map; and transforming the earth's field vector intobody coordinates to find errors in the “truth” source according to{tilde over (h)}_(i)=C_(L) ^(B)C_(N) ^(L{tilde over (h)}) ^(N), whereC_(L) ^(B) transforms a vector from the local-level frame (L-frame) tothe body frame (B-frame) and is the transpose of the attitude directioncosine matrix C, and C_(N) ^(L) accounts for the rotation in azimuth ofthe local-level frame with respect to north by the wander angle α, andis given by ${C_{N}^{L} = {\begin{bmatrix}{\cos\quad\alpha} & {\sin\quad\alpha} & 0 \\{{- \sin}\quad\alpha} & {\cos\quad\alpha} & 0 \\0 & 0 & 1\end{bmatrix}.}}\quad$
 31. A method according claim 30 furthercomprising taking partial differentials of {tilde over (h)}_(i)=C_(L)^(B)C_(N) ^(L){tilde over (h)}^(N) according to $\begin{matrix}{{\delta{\overset{\sim}{h}}_{i}} = {{\delta\quad C^{T}C_{N}^{L}{\overset{\sim}{h}}^{N}} + {C^{T}\delta\quad C_{N}^{L}{\overset{\sim}{h}}^{N}} + {C^{T}C_{N}^{L}\delta{\overset{\sim}{h}}^{N}}}} \\{{= {{\left\lbrack {C^{T}\left\{ \gamma \right\}} \right\rbrack C_{N}^{L}{\overset{\sim}{h}}^{N}} + {{C^{T}\left\lbrack {{- \left\{ ɛ \right\}}C_{N}^{L}} \right\rbrack}{\overset{\sim}{h}}^{N}} + {C^{T}C_{N}^{L}\delta{\overset{\sim}{h}}^{N}}}},}\end{matrix}\quad$ where γ is the attitude error vector which representsthe angular error of the L-frame relative to the B-frame, ε the angularposition error vector which represents the angular error of the L-framerelative to the earth frame (E-frame), and {V} represents theskew-symmetric matrix form of the enclosed vector v, defined by${\left\{ v \right\} \equiv {\begin{bmatrix}0 & {- v_{z}} & {\quad v_{y}} \\v_{z} & 0 & {- v_{x}} \\{- v_{y}} & v_{x} & 0\end{bmatrix}.}}\quad$
 32. A method according claim 31 furthercomprising: estimating corrections using a Kalman filter implementing a“psi-angle” inertial error model where attitude error states arecomponents of the angular error vector Ψ defined as Ψ=γ−ε, which resultsin δ{tilde over (h)}_(i)=C^(T){Ψ}C_(N) ^(L){tilde over(h)}^(N)+C^(T)C_(N) ^(L)δ{tilde over (h)}^(N), or by taking a reversedorder cross product, {Ψ}C_(N) ^(L){tilde over (h)}^(N)=Ψx (C_(N)^(L){tilde over (h)}^(N))=−(C_(N) ^(L){tilde over (h)}^(N))xΨ=−{C_(N)^(L){tilde over (h)}^(N)}Ψ, which results in δ{tilde over(h)}_(i)=−C^(T){C_(N) ^(L){tilde over (h)}^(N)}Ψ+C^(T)C_(N) ^(L)δ{tildeover (h)}^(N).
 33. A method according to claim 32 wherein determining adifference between a vector as measured by the magnetometer and thetruth reference vector comprises stating the differences asz_(k)=H_(k)x_(k)+v_(k), where z_(k) is a 3×1 vector of measurements attime t_(k), H_(k) is a 3×n measurement matrix, x_(k) is the statevector, and v_(k) is the measurement noise vector.
 34. A methodaccording to claim 33 further comprising: determining the measurementvector by subtracting δ{tilde over (h)}_(i)=−C^(T){C_(N) ^(L){tilde over(h)}^(N)}Ψ+C^(T)C_(N) ^(L)δ{tilde over (h)}^(N) from${{{\delta\overset{\sim}{h}} = {{\begin{bmatrix}{\Delta\quad{\overset{\sim}{h}}^{T}} & 0_{1 \times 3} & 0_{1 \times 3} \\0_{1 \times 3} & {\Delta\quad{\overset{\sim}{h}}^{T}} & 0 \\0_{1 \times 3} & 0_{1 \times 3} & {\Delta\quad{\overset{\sim}{h}}^{T}}\end{bmatrix}\quad\delta\quad I} - {L\quad\delta\quad{\overset{\sim}{h}}_{p}}}},}\quad$resulting in ${\begin{matrix}{z = {{\overset{\sim}{h} - {\overset{\sim}{h}}_{i}} = {{\delta\overset{\sim}{h}} - {\delta{\overset{\sim}{h}}_{i}} + v}}} \\{= {{\begin{bmatrix}{\Delta\quad{\overset{\sim}{h}}^{T}} & 0_{1 \times 3} & 0_{1 \times 3} \\0_{1 \times 3} & {\Delta\quad{\overset{\sim}{h}}^{T}} & 0 \\0_{1 \times 3} & 0_{1 \times 3} & {\Delta\quad{\overset{\sim}{h}}^{T}}\end{bmatrix}\quad\delta\quad I} - {L\quad\delta\quad{\overset{\sim}{h}}_{p}} + {C^{T}\left\{ {C_{N}^{L}\quad{\overset{\sim}{h}}^{N}} \right\}\psi} - {C^{T}C_{N}^{L}\delta{\overset{\sim}{h}}^{N}} + v}}\end{matrix},}\quad$ where v represents a vector of uncorrelatedmeasurement noise due to electromagnetic noise and random earth fieldmodeling errors such as quantization.
 35. A method according to claim 34further comprising generating a measurement mapping matrix H accordingto H=[H_(Ψ)0 . . . 0 H_(δl) H_({tilde over (h)}) _(p)H_({tilde over (h)}) _(N) ], where $\begin{matrix}{H_{\psi} = {C^{T}\left\{ {C_{N}^{L}\quad{\overset{\sim}{h}}^{N}} \right\}}} \\{{H_{\delta 1} = {\begin{bmatrix}{\Delta\quad{\overset{\sim}{h}}^{T}} & 0_{1 \times 3} & 0_{1 \times 3} \\0_{1 \times 3} & {\Delta\quad{\overset{\sim}{h}}^{T}} & 0 \\0_{1 \times 3} & 0_{1 \times 3} & {\Delta\quad{\overset{\sim}{h}}^{T}}\end{bmatrix}.}}\quad} \\{H_{{\overset{\sim}{h}}_{p}} = {- L}} \\{H_{{\overset{\sim}{h}}^{N}} = {{- C^{T}}C_{N}^{L}}}\end{matrix}$
 36. A magnetic compass compensation unit forcharacterizing distortions in the earth's magnetic field caused by avehicle, in which said unit is mounted, relative to an undisturbedmagnetic field of the earth, said unit comprising a processor configuredto receive measurements of the distorted magnetic field from amagnetometer, receive an orientation of the vehicle axes with respect tothe earth at a time corresponding to each magnetometer measurement,receive undistorted earth magnetic field data for the vicinity of thevehicle relative to the earth at the time corresponding to eachmagnetometer measurement, and characterize the distortions utilizing themagnetic field measurements, the orientations of the vehicle, and theundistorted earth magnetic field data.
 37. A magnetic compasscompensation unit according to claim 36, said processor receiving pitchand roll data from an attitude heading reference system.
 38. A magneticcompass compensation unit according to claim 36 wherein to characterizethe distortions, said processor calculates correction coefficients basedupon a magnetometer orientation with respect to the vehicle axes, themagnetic field as measured by the magnetometer, magnetic inclination anddeclination data, and vehicle orientation.
 39. A magnetic compasscompensation unit according to claim 38 wherein said processor isconfigured to calculate correction coefficients by estimating L_(e) andH_(pe) from the function H_(earth)=L_(e) * H_(meas)+H_(pe) where,H_(earth) is a three-dimensional vector representing the undisturbedEarth's magnetic field, H_(meas) represents a three-dimensional vectorfor the disturbed Earth's magnetic field as measured by magnetometer,L_(e) is a three by three (3×3) matrix of magnetic correctioncoefficients, and H_(pe) is a three-dimensional vector of magneticcorrection coefficients.
 40. A magnetic compass compensation unitaccording to claim 36 further configured to receive signalscorresponding to magnetic sensors oriented in at least two differentorientations.
 41. A magnetic compass compensation unit according toclaim 36 further configured to receive signals corresponding to threemagnetic sensors orthogonally oriented to one another.
 42. A magneticcompass compensation unit according to claim 36 further configured todetermine vehicle orientation based upon a received GPS tracking angle.43. A magnetic compass compensation unit according to claim 36 furtherconfigured to determine a vehicle orientation based upon a true headingcalculated utilizing GPS-augmented attitude and heading reference systemalgorithms.
 44. A magnetic compass compensation unit according to claim36 further configured to determine a vehicle orientation based upon atleast one of a true heading signal or a magnetic heading signal from aninertial reference system.
 45. A magnetic compass compensation unitaccording to claim 36 further configured to determine an orientation ofthe vehicle based upon a true heading received from a dual antenna GPS.46. A magnetic compass compensation unit according to claim 36 whereinto receive undistorted earth magnetic field data said processor accessesa table of magnetic inclination and declination data for a surface ofthe earth.
 47. A magnetic compass compensation unit according to claim46 wherein said processor is configured to transform undistorted earthmagnetic field data to correspond to axes of the magnetometer.
 48. Aprocessor programmed to generate a truth reference field vector, {tildeover (h)}_(i), from inertial data and a three dimensional map of theearth's magnetic field, determine a difference between a vector asmeasured by the magnetometer and the truth reference vector, and utilizethe difference to estimate corrections to magnetometer modelcoefficients.
 49. A processor according to claim 48 wherein saidprocessor comprises a Kalman filter for estimating corrections.
 50. Aprocessor according to claim 48 wherein said processor is programmed togenerate a truth reference field vector by relating the field measuredby a magnetometer {tilde over (h)}_(m) to the true earth's field {tildeover (h)} according to {tilde over (h)}_(m)=M{tilde over (h)}+{tildeover (h)}_(p), where M is a 3×3 magnetic permeability matrix and {tildeover (h)}_(p) is a 3×1 vector of field offset errors resulting frompermanent magnetization.
 51. A processor according to claim 50 whereinsaid processor is programmed to solve for the true earth's field, {tildeover (h)}, according to ${\begin{matrix}{\overset{\sim}{h} = {M^{- 1}\left( {{\overset{\sim}{h}}_{m} - {\overset{\sim}{h}}_{p}} \right)}} \\{= {L\left( {{\overset{\sim}{h}}_{m} - {\overset{\sim}{h}}_{p}} \right)}}\end{matrix},}\quad$ by assuming that the initial value of matrix L isequal to the identity matrix and the initial vector of field offseterrors, {tilde over (h)}_(p), is zero.
 52. A processor according toclaim 51 wherein said processor is programmed to: provide a linearizederror equation for the magnetometer-measured and compensated body-axisearth's field components according to $\begin{matrix}{{\delta\quad\overset{\sim}{h}} = {{\delta\quad{L\left( {{\overset{\sim}{h}}_{m} - {\overset{\sim}{h}}_{p}} \right)}} - {L\quad\delta\quad{\overset{\sim}{h}}_{p}}}} \\{{= {{\delta\quad L\quad\Delta\quad\overset{\sim}{h}} - {L\quad\delta\quad{\overset{\sim}{h}}_{p}}}},}\end{matrix}\quad$ where Δ{tilde over (h)} is defined as Δ{tilde over(h)}={tilde over (h)}_(m)−{tilde over (h)}_(p); and generate a separatemeasurement equation for each field component from the linearized errorequation according to δ{tilde over (h)}_(x)=Δ{tilde over(h)}^(T)δl_(r1)−l_(r1) ^(T)δ{tilde over (h)}_(p), δ{tilde over(h)}_(y)=Δ{tilde over (h)}^(T)δl_(r2)−l_(r2) ^(T)δ{tilde over (h)}_(p),and δ{tilde over (h)}_(z)=Δ{tilde over (h)}^(T)δl_(r3)−l_(r3)^(T)δ{tilde over (h)}_(p), where each equation is in the form of a rowvector multiplied by a column vector of magnetometer parameter errors toprovide scalars that are functions of the vectors.
 53. A processoraccording to claim 52 wherein said processor is programmed withmagnetometer errors according to${{\delta\quad\overset{\sim}{h}} = {{\begin{bmatrix}{\Delta\quad{\overset{\sim}{h}}^{T}} & 0_{1 \times 3} & 0_{1 \times 3} \\0_{1 \times 3} & {\Delta\quad{\overset{\sim}{h}}^{T}} & 0 \\0_{1 \times 3} & 0_{1 \times 3} & {\Delta\quad{\overset{\sim}{h}}^{T}}\end{bmatrix}\quad\delta\quad l} - {L\quad\delta\quad{{\overset{\sim}{h}}_{p}.}}}}\quad$54. A processor according to claim 48 wherein to generate a truthreference field vector, {tilde over (h)}_(i), from inertial data and athree dimensional map of the earth's magnetic field, said processor isprogrammed to: determine the earth's field vector in north/east/downframe components {tilde over (h)}^(N) based on one or more of a currentlatitude, longitude, altitude, and the map; and transform the earth'sfield vector into body coordinates to find errors in the “truth” sourceaccording to {tilde over (h)}_(i)=C_(L) ^(B)C_(N) ^(L){tilde over(h)}^(N), where C_(L) ^(B) transforms a vector from the local-levelframe (L-frame) to the body frame (B-frame) and is the transpose of theattitude direction cosine matrix C, and C_(N) ^(L) accounts for therotation in azimuth of the local-level frame with respect to north bythe wander angle α, and is given by ${C_{N}^{L} = {\begin{bmatrix}{\cos\quad\alpha} & {\sin\quad\alpha} & 0 \\{{- \sin}\quad\alpha} & {\cos\quad\alpha} & 0 \\0 & 0 & 1\end{bmatrix}.}}\quad$
 55. A processor according claim 54 wherein saidprocessor is programmed to take partial differentials of {tilde over(h)}_(i)=C_(L) ^(B)C_(N) ^(L){tilde over (h)}^(N) according to$\begin{matrix}{{\delta{\overset{\sim}{h}}_{i}} = {{\delta\quad C^{T}C_{N}^{L}{\overset{\sim}{h}}^{N}} + {C^{T}\delta\quad C_{N}^{L}{\overset{\sim}{h}}^{N}} + {C^{T}C_{N}^{L}\delta\quad{\overset{\sim}{h}}^{N}}}} \\{{= {{\left\lbrack {C^{T}\left\{ \gamma \right\}} \right\rbrack\quad C_{N}^{L}{\overset{\sim}{h}}^{N}} + {{C^{T}\left\lbrack {{- \left\{ ɛ \right\}}\quad C_{N}^{L}} \right\rbrack}\quad{\overset{\sim}{h}}^{N}} + {C^{T}C_{N}^{L}\delta\quad{\overset{\sim}{h}}^{N}}}},}\end{matrix}\quad$ where γ is the attitude error vector which representsthe angular error of the L-frame relative to the B-frame, ε the angularposition error vector which represents the angular error of the L-framerelative to the earth frame (E-frame), and {V} represents theskew-symmetric matrix form of the enclosed vector v, defined by${\left\{ v \right\} \equiv {\begin{bmatrix}0 & {- v_{z}} & v_{y} \\v_{z} & 0 & {- v_{x}} \\{- v_{y}} & v_{x} & 0\end{bmatrix}.}}\quad$
 56. A processor according claim 55 wherein saidprocessor is programmed to estimate corrections using a Kalman filterwhich implements a “psi-angle” inertial error model, the inertial errormodel having attitude error states that are components of an angularerror vector Ψ defined as Ψ=γ−ε, which results in an error model ofδ{tilde over (h)}_(i)=C^(T) {Ψ}C_(N) ^(L){tilde over (h)}^(N)+C^(T)C_(N)^(L)δ{tilde over (h)}^(N).
 57. A processor according to claim 56 whereinsaid processor is programmed to determine a difference between a vectoras measured by the magnetometer and a truth reference vector comprisesaccording to z_(k)=H_(k)x_(k)+v_(k), where z_(k) is a 3×1 vector ofmeasurements at time t_(k), H_(k) is a 3×n measurement matrix, x_(k) isthe state vector, and V_(k) is the measurement noise vector.
 58. Aprocessor according to claim 57 wherein said processor is programmed todetermine the measurement vector by subtracting δ{tilde over(h)}_(i)=−C^(T) {C_(N) ^(L){tilde over (h)}^(N)}Ψ+C^(T)C_(N) ^(L)δ{tildeover (h)}^(N) from ${{{\delta\overset{\sim}{h}} = {{\begin{bmatrix}{\Delta\quad{\overset{\sim}{h}}^{T}} & 0_{1 \times 3} & 0_{1 \times 3} \\0_{1 \times 3} & {\Delta\quad{\overset{\sim}{h}}^{T}} & 0 \\0_{1 \times 3} & 0_{1 \times 3} & {\Delta\quad{\overset{\sim}{h}}^{T}}\end{bmatrix}\quad\delta\quad l} - {L\quad\delta\quad{\overset{\sim}{h}}_{p}}}},}\quad$resulting in $\begin{matrix}{z = {\overset{\sim}{h} - {\overset{\sim}{h}}_{i}}} \\{= {{\delta\quad\overset{\sim}{h}} - {\delta{\overset{\sim}{h}}_{i}} + v}} \\{{= {{\begin{bmatrix}{\Delta\quad{\overset{\sim}{h}}^{T}} & 0_{1 \times 3} & 0_{1 \times 3} \\0_{1 \times 3} & {\Delta\quad{\overset{\sim}{h}}^{T}} & 0 \\0_{1 \times 3} & 0_{1 \times 3} & {\Delta\quad{\overset{\sim}{h}}^{T}}\end{bmatrix}\quad\delta\quad l} - {L\quad\delta\quad{\overset{\sim}{h}}_{p}} + {C^{T}\left\{ {C_{N}^{L}{\overset{\sim}{h}}^{N}} \right\}\quad\psi} - {C^{T}C_{N}^{L}\delta\quad{\overset{\sim}{h}}^{N}} + v}},}\end{matrix}\quad$ where v represents a vector of uncorrelatedmeasurement noise due to electromagnetic noise and random earth fieldmodeling errors such as quantization.
 59. A processor according to claim58 wherein said processor is programmed to generate a measurementmapping matrix H according to H=[H_(Ψ)0 . . . 0 H_(δl)H_({tilde over (h)}) _(p) H_({tilde over (h)}) _(N) ], where$\begin{matrix}{H_{\psi} = {C^{T}\left\{ {C_{N}^{L}{\overset{\sim}{h}}^{N}} \right\}}} \\{H_{\delta\quad l} = \begin{bmatrix}{\Delta\quad{\overset{\sim}{h}}^{T}} & 0_{1 \times 3} & 0_{1 \times 3} \\0_{1 \times 3} & {\Delta\quad{\overset{\sim}{h}}^{T}} & 0 \\0_{1 \times 3} & 0_{1 \times 3} & {\Delta\quad{\overset{\sim}{h}}^{T}}\end{bmatrix}} \\{H_{{\overset{\sim}{h}}_{p}} = {- L}} \\{H_{{\overset{\sim}{h}}^{N}} = {{- C^{T}}C_{N}^{L}}}\end{matrix}\quad$